Raw R code and data visualization for CASA0005

The main purpose for this assignment’s report is to find out the spatial distribution of COVID-19 new cases in London borough and to quantify the effect of ‘Lockdown’ policy on the suppress of the increase of new cases based on statistical analysis and spatial autocorrelation.

Data visualization on map

#import library 
library(sf)
library(tmap)
library(tmaptools)
library(tidyverse)
library(here)
library(janitor)
library(RColorBrewer)
library(sp)
library(spdep)

After completing loading libraries, loading all datasets needed for the analysis:

#read all dataset for the analysis
#London borough shapefile
Londonborough <-
  st_read(
    here::here(
      'data',
      'statistical-gis-boundaries-london',
      'ESRI',
      'London_Borough_Excluding_MHW.shp'
    )
  ) %>%
  st_transform(., 27700)
## Reading layer `London_Borough_Excluding_MHW' from data source `/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/statistical-gis-boundaries-london/ESRI/London_Borough_Excluding_MHW.shp' using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 7 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 503568.2 ymin: 155850.8 xmax: 561957.5 ymax: 200933.9
## projected CRS:  OSGB 1936 / British National Grid

Read csv file–>cumulative new cases in each time period, raw csv files are in Github repo

#cumulative new cases in 3.23-4.22 period
covid_323_422 <-
  read_csv(
    '/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-323-422.csv',
    na = c("NA", "n/a")
  ) %>%
  clean_names()
## 
## ─ Column specification ─────────────────────────────────────────────────
## cols(
##   area_code = col_character(),
##   new_cases = col_double(),
##   `population(10k)` = col_double(),
##   `new_cases per 10k population` = col_double()
## )
#cumulative new cases in 5.23-6.22 period
covid_523_622 <-
  read_csv(
    '/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases-523-622.csv',
    na = c("NA", "n/a")
  ) %>%
  clean_names()
## 
## ─ Column specification ─────────────────────────────────────────────────
## cols(
##   area_code = col_character(),
##   new_cases = col_double(),
##   population_10k = col_double(),
##   `new_cases_per_10k_523-622` = col_double()
## )
names(covid_323_422)
## [1] "area_code"                    "new_cases"                   
## [3] "population_10k"               "new_cases_per_10k_population"
names(covid_523_622)
## [1] "area_code"                 "new_cases"                
## [3] "population_10k"            "new_cases_per_10k_523_622"

Connect the CSV file data and shapefile

#join data with shapefile and new cases data in 3.23-4.22 period
Londonborough_merged <-
  Londonborough %>% left_join(covid_323_422, by = c('GSS_CODE' = 'area_code')) %>%
  distinct(GSS_CODE, NAME, new_cases,new_cases_per_10k_population) #choose which column to be used in the following analysis

#join data with shapefile and new cases data in 5.23-6.22 period
Londonborough_merged_523 <-
  Londonborough %>% left_join(covid_523_622, by = c('GSS_CODE' = 'area_code')) %>%
  distinct(GSS_CODE, NAME, new_cases, new_cases_per_10k_523_622) #choose which column to be used in the following analysis
#making maps
tmap_mode('view')
## tmap mode set to interactive viewing
breaks_323 = c(10, 15, 20, 25, 30, 35) 
tm1 <- tm_shape(Londonborough_merged) +
  tm_polygons(
    'new_cases_per_10k_population',
    breaks = breaks_323,
    alpha = 0.5,
    palette = 'PuBu',
    colorNA = 'white',
    title = 'New cases rate'
  ) + tm_layout(
    legend.position = c('left', 'bottom'),
    legend.outside = FALSE,
    legend.title.size = 1.2,
    legend.text.size = 0.75,
    frame = FALSE
  ) + tm_credits('(A) New cases rate in 3.23-4.22', position = c('left', 'top'), size = 1.2) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =c('right', 'bottom'), text.size = 0.6)
#view map 1
tm1
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
#saving map
tmap_save(tm1,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(323-422).png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#making map 2
breaks_523=c(0,0.5,1,1.5,2,2.5,3)
tm2 <- tm_shape(Londonborough_merged_523) +
  tm_polygons(
    'new_cases_per_10k_523_622',
    breaks = breaks_523,
    alpha = 0.5,
    palette = 'PuBu',
    colorNA = 'white',
    title = 'New cases rate'
  ) + tm_layout(
    legend.position = c('left', 'bottom'),
    legend.outside = FALSE,
    legend.title.size = 1.2,
    legend.text.size = 0.75,
    frame = FALSE
  ) +
  tm_compass(type = "arrow", position = c("right", "bottom")) + tm_credits('(B) New cases rate in 5.23-6.22',
                                                                           position = c('left', 'top'),
                                                                           size =
                                                                             1.2) +
  tm_scale_bar(position = c('right', 'bottom'), text.size = 0.6)
tm2
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm2,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/new_cases_rate(523-622).png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#combine two map
t=tmap_arrange(tm1,tm2,nrow=1)
tmap_save(t,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/New_cases_rate.png
## Resolution: 2100 by 2100 pixels
## Size: 7 by 7 inches (300 dpi)
t
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.

###Global Moran

#data cleaning for 3.23-4.22 data
london_exclude_city <-
  na.omit(Londonborough_merged) #exclude rows contain null data
#data cleaning for 5.23-6.22 data
london_exclude_city_523 <- na.omit(Londonborough_merged_523)

#creating polygon for 323 dataset and 523 dataset
neibour_323 <- poly2nb(london_exclude_city, queen = TRUE)
neibour_323[[1]]
## [1] 19 20 21 22
neibour_523 <- poly2nb(london_exclude_city_523,queen=TRUE)
neibour_523[[1]]
## [1] 19 20 21 22
#assign weight matrix for each neighbouring polygon using row standardisation
lw_323 <- nb2listw(neibour_323, style = "W", zero.policy = TRUE) #each neighbour is assigned a quarter of total weight
lw_523 <- nb2listw(neibour_523,style='W',zero.policy = TRUE)

#computing global Moran's I 
moran.test(london_exclude_city$new_cases_per_10k_population, lw_323)
## 
##  Moran I test under randomisation
## 
## data:  london_exclude_city$new_cases_per_10k_population  
## weights: lw_323    
## 
## Moran I statistic standard deviate = 3.1736, p-value = 0.0007528
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.35308476       -0.03225806        0.01474307
moran.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
## 
##  Moran I test under randomisation
## 
## data:  london_exclude_city_523$new_cases_per_10k_523_622  
## weights: lw_523    
## 
## Moran I statistic standard deviate = 2.5898, p-value = 0.004802
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.28295263       -0.03225806        0.01481388

###Local Moran visualization

#local moran value for data 3.23-4.22
local_moran_323 <- localmoran(london_exclude_city$new_cases_per_10k_population, lw_323)
summary(local_moran_323)
##        Ii                 E.Ii              Var.Ii            Z.Ii         
##  Min.   :-0.419544   Min.   :-0.03226   Min.   :0.1102   Min.   :-0.71017  
##  1st Qu.:-0.002326   1st Qu.:-0.03226   1st Qu.:0.1664   1st Qu.: 0.08039  
##  Median : 0.149082   Median :-0.03226   Median :0.2155   Median : 0.36073  
##  Mean   : 0.353085   Mean   :-0.03226   Mean   :0.2460   Mean   : 0.82412  
##  3rd Qu.: 0.619470   3rd Qu.:-0.03226   3rd Qu.:0.2974   3rd Qu.: 1.59790  
##  Max.   : 2.209853   Max.   :-0.03226   Max.   :0.4612   Max.   : 4.11139  
##    Pr(z > 0)        
##  Min.   :0.0000197  
##  1st Qu.:0.0550665  
##  Median :0.3595183  
##  Mean   :0.3109317  
##  3rd Qu.:0.4679710  
##  Max.   :0.7612010
#local moran value for data 5.23-6.22
local_moran_523 <- localmoran(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
summary(local_moran_523)
##        Ii                E.Ii              Var.Ii            Z.Ii         
##  Min.   :-0.40785   Min.   :-0.03226   Min.   :0.1106   Min.   :-0.80748  
##  1st Qu.:-0.05669   1st Qu.:-0.03226   1st Qu.:0.1670   1st Qu.:-0.05456  
##  Median : 0.16490   Median :-0.03226   Median :0.2164   Median : 0.43173  
##  Mean   : 0.28295   Mean   :-0.03226   Mean   :0.2470   Mean   : 0.68948  
##  3rd Qu.: 0.44756   3rd Qu.:-0.03226   3rd Qu.:0.2986   3rd Qu.: 1.17418  
##  Max.   : 1.39046   Max.   :-0.03226   Max.   :0.4632   Max.   : 3.48158  
##    Pr(z > 0)        
##  Min.   :0.0002492  
##  1st Qu.:0.1213157  
##  Median :0.3343064  
##  Mean   :0.3309740  
##  3rd Qu.:0.5217552  
##  Max.   :0.7903052
#plot local moran map
#There are 5 columns of data.
#copy some of the columns (the I score (column 1) and the z-score standard deviation (column 4))
#the z-score (standardised value relating to whether high values or low values are clustering together)
#change local_moran type
local_moran_tibble_323 <- as_tibble(local_moran_323) #change to dataframe
local_moran_tibble_523 <- as_tibble(local_moran_523)

london_exclude_city <-
  london_exclude_city %>% mutate(local_moran_I = as.numeric(local_moran_tibble_323$Ii)) %>% mutate(local_moran_z =as.numeric(local_moran_tibble_323$Z.Ii))

london_exclude_city_523 <-
  london_exclude_city_523 %>% mutate(local_moran_I = as.numeric(local_moran_tibble_523$Ii)) %>% mutate(local_moran_z =
                                                                                                         as.numeric(local_moran_tibble_523$Z.Ii))

#setting color
MoranColours <- rev(brewer.pal(8, "RdBu"))

#plot a map (in Rmd documend, tmap mode changes to viewing, but raw code is plotting mode)
tmap_mode('view')
## tmap mode set to interactive viewing
#plotting local moran map for 3.23-4.22
tm_323 <- tm_shape(london_exclude_city) +
  tm_polygons(
    'local_moran_I',
    alpha = 0.5,
    palette = MoranColours,
    title = 'Local Moran I',midpoint=NA
  ) + tm_layout(
    legend.position = c('left', 'bottom'),
    legend.outside = FALSE,
    legend.title.size = 1.2,
    legend.text.size = 0.75,
    frame = FALSE
  ) + tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
                                                                                  c('right', 'bottom'), text.size = 0.6)+tm_credits('(A) Local Moran in 3.23-4.22',position = c('left','top'),size = 1.2)

tm_323
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm_323,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_323.png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#plotting local moran map for 5.23-6.22 dataset
tm_523 <- tm_shape(london_exclude_city_523) +
  tm_polygons(
    'local_moran_I',
    alpha = 0.5,
    palette = MoranColours,
    title = 'Local Moran I',midpoint=NA
  ) + tm_layout(
    legend.position = c('left', 'bottom'),
    legend.outside = FALSE,
    legend.title.size = 1.2,
    legend.text.size = 0.75,
    frame = FALSE
  ) +tm_compass(type = "arrow", position = c("right", "bottom")) +tm_scale_bar(position =
                                                                                 c('right', 'bottom'), text.size = 0.6)+tm_credits('(B) Local Moran in 5.23-6.22',position = c('left','top'),size = 1.2)

tm_523
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(tm_523,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_523.png
## Resolution: 2389.896 by 1845.268 pixels
## Size: 7.966321 by 6.150895 inches (300 dpi)
#combine two map
t_c=tmap_arrange(tm_323,tm_523,nrow=1)
t_c
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
## Credits not supported in view mode.
## Compass not supported in view mode.
## legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
tmap_save(t_c,'/Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png')
## Map saved to /Users/ziyichen/Desktop/SDSV/GIS/Data_for_CW/Code&data/data/local_moran_combined.png
## Resolution: 2100 by 2100 pixels
## Size: 7 by 7 inches (300 dpi)

###Geary’s test result

#geary's test
#geary's test for 3.23-4.22 data
geary.test(london_exclude_city$new_cases_per_10k_population, lw_323)
## 
##  Geary C test under randomisation
## 
## data:  london_exclude_city$new_cases_per_10k_population 
## weights: lw_323 
## 
## Geary C statistic standard deviate = 2.9614, p-value = 0.001531
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.63209605        1.00000000        0.01543425
#geary's test for 5.23-6.22 data
geary.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
## 
##  Geary C test under randomisation
## 
## data:  london_exclude_city_523$new_cases_per_10k_523_622 
## weights: lw_523 
## 
## Geary C statistic standard deviate = 2.1708, p-value = 0.01497
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##        0.73102781        1.00000000        0.01535192

###Getis Ord General G result

#Getis Ord General G for 3.23-4.22
globalG.test(london_exclude_city$new_cases_per_10k_population,lw_323)
## Warning in globalG.test(london_exclude_city$new_cases_per_10k_population, :
## Binary weights recommended (sepecially for distance bands)
## 
##  Getis-Ord global G statistic
## 
## data:  london_exclude_city$new_cases_per_10k_population 
## weights: lw_323 
## 
## standard deviate = 2.5827, p-value = 0.004901
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       3.347784e-02       3.225806e-02       2.230541e-07
#Getis Ord General G for 5.23-6.22
globalG.test(london_exclude_city_523$new_cases_per_10k_523_622,lw_523)
## Warning in globalG.test(london_exclude_city_523$new_cases_per_10k_523_622, :
## Binary weights recommended (sepecially for distance bands)
## 
##  Getis-Ord global G statistic
## 
## data:  london_exclude_city_523$new_cases_per_10k_523_622 
## weights: lw_523 
## 
## standard deviate = 1.525, p-value = 0.06362
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       3.336104e-02       3.225806e-02       5.230748e-07

###Descriptive statistics summary

library(ggplot2)
#histrogram
#only use rate of new cases to plot histrogram
covid_323_hist <- hist(covid_323_422$new_cases_per_10k_population,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 3.23-4.22',angle=45,ylim=c(0,15))

mean_323 <- mean(covid_323_422$new_cases_per_10k_population)
var_323 <- var(covid_323_422$new_cases_per_10k_population)
std_323 <- sd(covid_323_422$new_cases_per_10k_population)
mean_323
## [1] 23.5875
var_323
## [1] 25.43532
std_323
## [1] 5.043344
covid_523_hist <- hist(covid_523_622$new_cases_per_10k_523_622,col='Dark blue',density=10,xlab='Rate of new cases',main='Rate of new cases during 5.23-6.22',angle=45,ylim=c(0,15))

mean_523 <- mean(covid_523_622$new_cases_per_10k_523_622)
var_523 <- var(covid_523_622$new_cases_per_10k_523_622)
std_523 <- sd(covid_523_622$new_cases_per_10k_523_622)
mean_523
## [1] 1.771875
var_523
## [1] 0.2949899
std_523
## [1] 0.5431297